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Abstract 

Finite  difference  methods  for  plate  buckling  were 
investigated  for  various  two-dimensional  plates.  The  plates 
were  clamped  and  simply  supported.  The  rate  of  convergence 
to  buckling  coefficients,  considering  virtual  work  and  equi- 
librium, have  been  compared.  A trigonometric  function  (based 
on  long  plate  theory)  was  incorporated  into  the  finite 
difference  approximations  of  the  virtual  work  expression. 
Additional  trigonometric  parameters  were  also  considered 
in  the  virtual  work  equation.  Little  difference  in  accuracy 
was  found  between  results  obtained  from  conventional  (poly- 
nomial) virtual  work  versus  the  equilibrium  approach.  Notice- 
able improvement  is  obtained  by  using  a trigonometric  function 
(long  plate  assumption)  in  the  virtual  work  expression.  If 
one  varies  the  trigonometric  function  (not  based  on  long 
plate  theory) , a band  of  results  is  created  from  use  of  the 
virtual  work  equation.  Inaccuracies  occurred  at  low  degrees 
of  freedom,  for  clamped  plates  at  all  length  to  width  ratios, 
due  to  insufficient  representation  of  the  boundaries.  Two 
dominant  reasons  for  poor  boundary  representation  seem  to 
be?  large  mesh  size  at  low  degrees  of  freedom  (mesh  sizes 
were  held  constant)  and  the  trigonometric  function  (poly- 
nomial for  conventional  finite  difference)  failing  to  satisfy 
the  kinematic  boundary  conditions. 
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COMPARING  TRIGONOMETRIC  AND  CONVENTIONAL 


FINITE  DIFFERENCE  APPROXIMATIONS  FOR 
PLATE  BUCKLING 


I.  Introduction 


Background 

The  plate  is  one  of  the  most  effective  aerospace  struc- 
tural elements.  Many  plate  characteristics  are  associated 
with  thinness,  thus  thinness  becomes  an  important  criteria 
for  optimizing  the  weight  considerations  of  the  overall 
structure.  The  concern  for  optimization  involves  stress 
calculations  which  consider  equilibrium.  Such  calculations 
are  not  in  themselves  the  end  results,  for  an  analyst  must 
go  one  step  further  and  determine  the  buckling  parameters  of 
the  plate  structure.  This  thesis  deals  with  a numerical 
approach  to  plate  buckling  using  a recently  developed  finite 
difference  technique  [1] . 

The  theoretical  buckling  stress  of  a flat  plate  is  the 
stress  at  which  an  exchange  of  stable  equilibrium  occurs 
between  the  straight  and  slightly  bent  form  of  the  plate. 

It  marks  the  region  where  continued  load  application  accel- 
erates the  growth  of  deflections  perpendicular  to  the  plane 
of  the  plate.  The  importance  of  this  stability  condition 
lies  in  the  fact  that  buckling  initiates  the  processes  lead- 
ing to  plate  failure  [2] . Because  of  the  necessity  to 
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incorporate  difficult  boundary  conditions  [3] , very  few 
plate  buckling  problems  of  a practical  engineering  nature 
have  been  solved  by  the  so-called  classical  techniques 
(series  solution) . Apparently  the  earliest  solution  to  a 
flat  plate  stability  problem  was  presented  by  Bryan  [4]  in 
1891  [5] . Since  then,  many  additional  solution  techniques 
have  been  developed  including  Navier  Solutions  [6] , use  of 
Fourier  Series  in  the  energy  equation  [7] , Galerkins  Method, 
and  other  approaches  [8-14].  All  are  used  to  derive  approx- 
imate answers,  but  they  are  lengthy  and  tedious  techniques. 
The  advent  of  the  high-speed  computer  and  subsequent  rapid 
development  of  numerical  methods  has  allowed  the  engineer 
the  opportunity  to  address  more  precise  conceptual  represen- 
tations of  actual  physical  structures. 

The  succeeding  sections  of  this  thesis  are  directed 
towards  the  development  of  numerical  approaches  to  buckling 
using  a CDC6600  computer.  The  most  efficient  numerical 
technique,  referred  to  as  rapid  convergence,  is  defined 
herein  as  that  technique  which  provided  useful  answers 
incorporating  the  fewest  degrees  of  freedom.  Rapid  con- 
vergence is  a definite  concern  for  buckling  considerstions . 

The  evaluation  of  the  convergence  rate  developed  from 
finite  difference  approximations  in  which  the  buckling  mode 
shape  characteristics  are  incorporated  into  analogous  grid 
spacing  is  of  prime  interest.  These  characteristics  were 
presented  by  M.  Stein  and  J.  Housner  in  1975  [1]  when  they 
introduced  a trigonometric  finite  difference  procedure.  A 
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comparison  between  the  new  trigonometric  technique  and  the 
finite  difference  approximation  using  the  common  equilibrium 
equation  is  presented  in  this  thesis. 


Purpose 

The  purpose  of  this  thesis  is  to  compare  buckling 
coefficients  determined  by  trigonometric  and  conventional 
finite  difference  methods  for  two-dimensional  plates,  incor- 
porating both  virtual  work  and  equilibrium  equations.  The 
development  of  rapid  convergence  is  of  primary  interest. 


General  Procedure 

In  the  approach  followed,  the  virtual  work  equation 
was  evaluated  using  finite  difference  approximations  of 
partial  derivatives.  This  required  the  separation  of  a 
plate's  domair  into  given  sets  of  nodes  which  were  then 
used  in  determining  each  partial  derivative  expression. 
Results  were  then  placed  into  the  virtual  work  equation. 
Integration  was  performed  leading  to  the  solution  of  a set 
of  algebraic  expressions.  Simple  support  and  clamped 
boundary  conditions  were  investigated.  Each  plate's  aspect 
ratio  (length/width=A/B)  was  varied  from  0.5  to  5.  Previous 
work  [1]  presented  the  results  for  trigonometric  finite 
difference  considering  square  plates  and  very  long  plates 
(A/B=5) . This  thesis  considered  the  borderline  between 
plate  and  slab  theory  by  investigating  aspect  ratios 
of  2 and  3.  Solutions  for  both  trigonometric  and  conven- 
tional finite  difference  approximations  of  the  virtual 


work  equation  have  been  obtained  and  compared  to  the  conven- 
tional finite  difference  approximations  of  the  equilibrium 
equation. 


! 


1 


II.  Theory 


Assumptions 

The  buckling  analysis  of  rectangular  plates  has  been 
carried  out  with  the  following  assumptions. 

1.  Linear  Theory  of  Elasticity  [6,15,16] 

a)  Plane  sections  remain  plane. 

b)  All  displacements  are  related  to  the  middle 
surface . 

c)  The  plate  is  perfectly  flat  prior  to  buckling. 

2.  Isotropic  material 

3.  Homogeneous  material 


4.  In-plane  load  Nx=constant  and  Nxy=N^=0 


5.  The  nonlinear  strain-displacement  relationships 


used  to  obtain  the  virtual  work  equation  are 


1 , ,2 

e =u  +-=-(w  ) 

x ,x  2'  ,x 


e =v  +i(w  ) ^ 

y *y  2V  ,y 


(2-1) 


e =u  +v  +w  w 
xy  ,y  ,x  ,x  ,y 


where  e ,e  , and  e are  the  strains  and  u,  v,  and  w are  the 
x'  y xy 


displacements  in  the  x-,y-,  and  z-  directions,  respectively 
[ 15]  . The  virtual  work  equation  is  developed  in  Appendix  A. 

6.  Thickness  of  plate  is  constant  and  the  plate  is  thin. 


Governing  Equations 

Equilibrium  Differential  Equation.  The  general  form 
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of  the  differential  equation  describing  the  slightly  bent 
equilibrium  configuration  of  an  initially  flat  plate  was 
derived  by  Stowell  [17]  and  reduced  to  equation  (2-2)  for 


the  elastic  case  in  reference  [ 2 ] . 

4 t 

V w = - —(a  w +2a  w +a  w ) 
D x ,xx  xy  ,xy  y ,yy 


(2-2) 


-to  = N 

X X 


to  = N =0 
xy  xy 


t0y  - Hy  « 0 (2-3) 


The  force  in  the  plane  of  the  plate  is  assumed  constant 
(assumption  No.  4)  and  is  thus  not  dependent  on  the  deflec- 
tion nor  is  it  a function  of  displacement.  This  makes 
equation  (2-2)  linear.  Hence  it  becomes 

(2-4) 

where  D is  the  flexural  rigidity.  Equation  (2-5)  defines 


DV  w = N w 

X , XX 


D for  an  isotropic  homogeneous  material. 

D = Et3/12 (1-u2) 


(2-5) 


Virtual  Work  Equation.  The  virtual  work  of 
the  plate  during  buckling  may  be  expressed  as 

6U  = 6V 

where 

fb  fa 


6U 


Uo(MxSu'xx+V^+2,Vw,xy,axdy 


(2-6) 


(2-7) 


6V  = [ [ (N  w 6w  ) dydx 
JoJo  x ' ' 


(2-8) 


The  sign  conventions  of  the  bending  moments  are  given  in 
figure  2-1.  Equation  (2-6)  is  developed  in  Appendix  A. 

6 


- 


III.  Numerical  Techniques 


In  the  nineteenth  century,  Boole  and  others  formalized 
the  finite  difference  method  for  solving  differential 
equations  [5] , thus  providing  the  means  to  reduce  the  con- 
tinuum to  a system  with  finite  degrees  of  freedom.  Mario 
G.  Salvadori  in  1949  [18]  presented  a comprehensive  paper 
concerning  the  difference  method  and  applicable  extrapola- 
tion procedures.  Presently,  many  excellent  texts  [19-25] 
may  be  found  covering  finite  differences.  This  thesis  deals 
with  the  central  difference  approximation.  Most  texts 
develop  full-station  central  difference,  however,  very 
little  is  written  on  the  use  of  half-stations.  This  section 
presents  a comparison  between  half-  and  full-station  central 
difference,  while  also  indicating  the  basic  steps  for 
coupling  the  finite  difference  method  with  plate  buckling. 

A new  trigonometric  approach  is  compared  to  the  standard 
polynomial  central  difference  approximation,  and  a plate 
grid  system  is  developed  to  show  how  boundary  conditions 
have  been  incorporated. 

Basic  Approach 

Plate  buckling,  as  considered  in  this  thesis,  makes 
use  of  the  central  finite  difference  method.  The  following 
steps  indicate  how  the  method  was  basically  employed  in 
the  analysis. 

1.  Discrete  grid  points  were  selected  to  approximate 
the  governing  partial  differential  equation  within 
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the  plate. 

2.  The  central  finite  difference  approximation  was 
applied  at  each  grid  and/or  reference  point. 

3.  Each  partial  differential  equation  was  reduced  to 
linear  algebraic  form. 

4.  The  resulting  equations  were  then  solved  and 
buckling  coefficient  determined. 

Half-Station  versus  Full-Station 

For  the  derivation  of  finite  difference  expressions,  it 
is  necessary  to  work  with  grid  arrangements.  Finite  dif- 
ference grids  for  a one-dimensional  case  are  shown  in  Figure 
3.1  and  Figure  3.2.  Grid  points  (indicated  by  circles) 
represent  discretely  defined  functional  values  while  ref- 
erence points,  indicated  by  an  x in  Figure  3.2,  locate  the 
position  of  the  functional  value  desired.  Grid  points  can 
be  coincident  with  reference  points  but  this  is  not  manda- 
tory. The  degrees  of  freedom  of  the  system  are  the  specific 
values  of  the  unknown  functions  located  at  grid  points.  The 
spacing  between  mesh  lines  was  held  constant.  The  grid 
lines  in  x-  and  y-  directions  are  represented  by  the  symbols 
hx  and  hy  respectively. 

To  develop  half-station  central  difference  equations, 
in  which  the  reference  point  is  located  between  grid  points, 
one  can  evaluate  a function  g,  at  point  (Fig  3.2)  using 
a Taylor  series  expansion.  This  of  course,  implies  continuity 
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Figure  3.1  Finite  Difference  Grid 


mesh 

line 


Rj  +J 


i ; 

. h L h ' 
| 2 T’  2 ; 
1 I 


Figure  3.2  Half-station  Finite  Difference  Grid 


conditions  for  each  functional  relationship.  The  basic 
Taylor  series  can  be  stated  as 

9(at-3/2)  “*9+ctii9,+3T  9ij9' '+TT  ai j9' ' '+*  * * * (at  Rj)  I3'1’ 
where  primes  denote  order  of  derivatives  and 

°ij  = (xat  -3/2  - xat  Rj>  <3-2> 

and  = h (a  constant  mesh  size) . It  is  now  possible  to 
incorporate  the  above  expressions  for  specific  locations. 


„ 3h,.l  f3hl2„,,  1 f3hl3  ,,,.1  (3hl 

g-3/2-*~T^  +2  [j-J  9 -r  In  9 +irlr-J 


4 iv  ,, 

g -...  (3-3) 
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g_ 


1/2 


h . 1 (hi  2 ,,  1 

=g-og'+-l-l  g -- 


•4 


24 


'*'4  iv 

g . 


g+l/2=g+2g'+l 


V2 


g',+M!)3g'"+lr 


43/2 


. 3h  ljLl 

=g+T-^'+o 


3h 


'g"4 


3h 


•4 


24 


3h 


iv, 

g +. . . 


(3-4) 


(3-5) 


(3-6) 


Equation  (3-4)  is  subtracted  from  equation  (3-5)  yielding 


i Yi 

g'=h  (g+l/2  ' g-l/2}  " 2l  g 


I I • 


(3-7) 


where  the  half-station  representation  of  the  first  derivative 
is 

'=-  ~ g_i  /o)  (3-8) 


g'=h  (g+l/2 


•1/2 ; 


h 

24 


and  Ur  g'1'  represents  the  error  approximation  [5],  On  the 


other  hand,  if  equations  (3-4)  and  (3-3)  are  added  to  equa- 
tions (3-5)  and  3-6) , respectively,  one  obtains  the  follow- 
ing expressions 


g+l/2+g-l/2  2g+ 


12 


iv 


(3-9) 


9+3/2+g-3/2  2g+ l 2 


3h 


^2 


1 ( 3h 


12  2 


.iv 


(3-10) 


One  can  solve  equation  (3-9)  for  2g  and  substitute  into 
equation  (3-10) 


g+3/2“g+l/2_g-l/2+g-3/2“2h2g' '+12  glV 


(3-11) 


This  last  equation  yields  the  second  derivative  based  upon 
half-stations 


. , _ 1 , „ . . 5h  iv 

9 ^T(g+3/2'g+l/2~g-l/2  g-3/2)-24  9 


(3-12) 
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Thus,  using  half-station  central  difference  for  the  first 
and  second  derivatives  and  removing  higher  order  terms,  one 


arrives  at 


9'  h (g+l/2_g-l/2) 


(3-13) 


g"  2h2  (g+3/2~g+l/2_g-l/2+g-3/2) 


(3-14) 


Full-station  central  difference  developed  in  Ref  [5]  is 
stated  here  for  completeness. 


2 

n iv 


g"  h2  (g+l_2g0+g-l)  "12  g 


(3-15) 


(3-16) 


It  becomes  apparent  that  first  order  derivatives  are  best 

approximated  between  grid  points  (half-station)  comparing 
h2  h2 

bound  error  g'*'  versus  g—  g'1'  yet,  second  derivatives 
can  be  more  accurately  determined  at  the  grid  points  (full- 
station)  if  the  higher  order  terms  become  a consideration 


h2  iv 


5h  _iv | 


12  g versus  2 T g 


Trigonometric  versus  Polynomial 


Therefore,  one  can  let  g(x)=P  where 

P = Ax2  + Bx  + C (3-17) 


A,  B,  and  C are  constants.  Now,  evaluate  (3-17)  at  points 
h,  o,  and  -h,  respectively. 


g(h)  = Ah2  + Bh  + C = g+1 

(3-18) 

g(o)  = C = gQ 

(3-19) 

g(-h)  = Ah2  - Bh  + C = g 

(3-20) 

where  g+1  etc.  are 

in  terms  of  full  station  central 

differ- 

ence  notation. 

Add  equations 

(3-18)  and  (3-20) 

g+1  + g_x  = 2 Ah 2 + 2C 

(3-21) 

Next  subtract 

2 times  equation  (3-19) 

2Ah2  = g+1  - 2go  + g_x 

(3-22) 

and  rearrange 

2A  = [g+l  " 2go  + g-l] 

(3-23) 
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where  the  right-hand  side  of  equation  (3-23)  is  the  defini- 
tion of  the  second-order  full  station  central  difference 


2 

(3-16) . One  can  now  take  the  second  derivative  (D  ) of 
equation  (3-17)  with  respect  to  x. 

D2g (x)  = 2A  (3-24) 

Thus,  equation  (3-24)  equals  equation  (3-23)  showing 

D2g(x)  = jT  [g+1  - 2go  + g_x]  (3-25) 

Consequently,  for  this  second  order  system  the  polynomial 
approximation  of  equation  (3-17)  becomes  an  equality. 
Conventionally  the  polynomial  approximation  was  used  with 
the  central  difference  technique.  However,  in  1975  M.  Stein 
and  J.  Housner  [1]  decided  to  let  the  function  g(x)  be  deter- 
mined by  a trigonometric  expression  (since  the  buckling  mode 
shape  is  normally  trigonometric  in  nature)  and  they  developed 
a trigonometric  finite-difference  solution  to  gain  a faster 
convergence  rate.  Equation  (3-26)  allows  g(x)  to  be  approxi- 
mated by  a trigonometric  expression. 

TT  (X— X ) TT(X— X ) 

g(x)  = + T2  sin  ~x + T3  COS  — X (3-26) 

The  symbol  represents  a wavelength  parameter  [1] . 

Earlier  paragraphs  showed  the  first  derivative  could  be 
found  by  half-station  central  difference,  as 

Dg(x)  = ~ (9+1/2  “ 9_i/2 ) (3-27) 

Through  the  use  of  equation  (3-26)  the  formation  of  an 
equality  can  be  determined  as  follows 


where 


(3-29 


2\  sin 
x 


Equation  (3-26)  is  repeated  for  convenience  and  its  develop 
ment  is  indicated  in  a step-by-step  fashion. 


tt  (x-x  ) 


tt(x-x  ) 


Take  derivative  of  equation  (3-26)  with  respect  to  x 


it  (x-x  ) 


it  (x-x_) 


Evaluate  equation  (3-30)  at  x=x 


where 


Dg(xQ)  = T2  tt/Xx 


T2 = y)g(xo)/lT 


Now  evaluate  equation  (3-26)  at  x 


n = Ti  + T2  sin 


f-  [ 

x 


X + ^ - X 
O 2 


+ T3  cos  b [xo  + £ - xo] 
/2  = Ti  + T2  sin  +T3  cos  2^- 
/2  = T1  + T2  Sln  [j  + T3  C 


E:  cos(-u)  = cos  u and  sin(-u) 


5.  Equation  (3-32)  is  substituted, 


^ Ax  Trh 

g+l/2  - g-l/2  = IT  Dg(xo>  sin  2T 


6.  Rearrange  equation  (3-37). 


(3-37) 


Dg(xQ)  = 


_ , irh 

2ix  Sln  2T 


-1/2  " g-l/2 


(3-38) 


Dg  (x 


o’  = fi  [9- 


1/2  ” g-l/2 


( 3-38a) 


The  reader  can  observe  in  equation  (3-38)  that  as  the  wave- 
length A^  approaches  infinity,  h approached  h,  since  the 
sin  (irh/2A  ) approached  (irh/2A  ) , thus  reducing  the  trigo- 
metric  expression  to  the  conventional  difference  expression. 
Both  trigonometric  and  polynomial  finite  difference  are 
used  to  approximate  the  virtual  work  equation. 


Mesh  Arrangement  in  the  Virtual  Work  Equation 

In  order  to  acquaint  the  reader  more  fully  with  the  use 
of  half-station  for  first  order  partial  derivatives  and  full- 
station  for  second  order  partials  an  example  follows, 
employing  the  trigonometric  technique.  The  virtual  work 
equation  (2-6)  is  repeated  here  for  completeness. 


LL  (MxSw,xx+MySw,yy+2HxySw.xy)dXdy=LL<Nx”.xSw,x)d''dx  {2~6) 
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Partial  derivatives  are  replaced  by  the  following  trigo- 
metric  central  difference  expressions 

<",xx>i,j  * <wi+l  . - 2wij  - ui-l  >/Sx 

» 1 fJ 


(w  ) . . = (w.  ...  - 2w. . + w.  . , )/n 

,yy  i,J  i»l+l  il  i,3-l  i 


(3-39) 


(w,xy>ij  = ' wi,j+l  ' Wi+l.j  + wij)/hxhy 

A A 

where  h and  h have  been  defined  in  equation  (3-29) . 
x y 

Equation  (2-6)  can  now  be  stated  as 

6u  - h*\  ?j=1  L [5xi5y:j[Dll<Wi+l"i  - 2”i: 


+ wi-l,j)/hx  + D12(”i,j+1  - 2uij  +Ui,:Kl)/ 
Kxfiy]  [5wi+l,j  - 2«“ij  * 5wi-l.j] 

+ SS  [D12<Wi«.J  ' 2”iJ  + 

+ D22(wi,j+1  - 2“ij  + wi/j-l*^y]  [6”i,j+l 
- 2Swi:j  + 6wljj_1]  + 2Vny.  [2D66(ui+l,j+l 


- wi,j+l  ‘ wi+l,j  + wij,<5wi+l,j+l  ‘ Swi,j+1 


- 5wi+l,j  + 


(3-40) 


where  integers  N and  M are  the  total  number  of  finite 
difference  stations  in  the  x-  and  y-  directions  and  the 
£ and  n coefficients  equal 
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Cxi  = o when  i < 3^  or  i > I 


= 1/2  when  i = ^ or  i = 1^ 


(3-41) 


= 1 when  1^  < i < i 


Cxj  = 0 when  j < JA  or  j > J, 


1/2  when  j = J4  or  j = J2 


= 1 when 


J4  < j < J2 


(3-42) 


nxi  = 0 when  1 < Ii  or  1 1 I3 


= 1 when  I1  £ i < I3 


nyj  = ° when  j < or  j > J 


- u2 


= 1 when 


(3-43) 


(3-44) 


J4  i 3 ' J2 

I1  and  I3  designate  the  grid  row  of  boundaries  1 and  3 while 
J2  and  are  column  designators  for  boundaries  2 and  4 
(Fig*  3.4) . The  right  hand  side  of  equation  (2-6)  can  also 
be  expanded  incorporating 

(3-45) 


w,x  - (wi+l,j  - "ij'/hx 

w,y  = <ui,j+i  - wij)/fiy 


(3-46) 


W>y  Sw,x  (wi,j+l  ■ wi j + ”i+l,j+l  ■ "i+1, j* (Swi+l, j 


- 6w . . + 6w. 


i+l, 3*1  - Swi,j+l,/4hxhy  <J-47> 
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and  a similar  arrangement  for  w Sw  . To  fully  appreciate 

/ x , y 

the  development  of  the  above  equations,  the  reader's  atten- 
tion is  directed  to  grid  point  7 of  a twenty-five  degree 
grid  arrangement  in  Figure  3.4.  The  numerical  approximation 
w and  its  variation  6w  are  evaluated  at  reference  point 

9 X fX 

A using  half-station  central  difference.  Thus  grid  point 
12  is  incorporated  with  grid  point  7. 

\x  {u,x  ■ (w12  - V(Sw12  • Su7)/fix  (3-43) 
The  same  procedure  is  applied  to  w 6w  in  the  y-direction. 
A problem  occurs  when  two  different  directional  derivatives 
(i.e.  w 6w  ) become  involved.  This  is  due  to  the  fact 

/ y »x 

that  derivative  functions  must  be  incorporated  into  the 

volume  integral  using  constant  quantities  (h  h (£  £ )). 

x y xi  yj 

In  the  case  referred  to  above  separate  ordinate  values 
(L  and  M)  would  be  obtained  at  reference  points  A and  C 
(Figure  3.5),  respectively 


■0  VetWME 

- — ~ ^ 

Vo* .om  i 

. ^ . 



Figure  3.5  Half-Station  Differences  of  w and  6w 

»y  rx 


Consequently,  two  ordinate  values  are  obtained,  but  they  are 
associated  with  two  distinct  volumes  for  the  same  functional 


w 6w  . In  order  to  maintain  a half-station  technique 
» x , y 

and  to  associate  the  product  of  two  separate  directional 
derivatives  with  a particular  volume,  a reference  point  is 
located  at  position  E and  is  evaluated  in  the  following 
manner. 


Figure  3.6  The  Position  of  a Reference  Point  to 
Evaluate  a Mixed  Derivative. 


The  half-station  central  difference  of  point  E is 

w = (wB  - «A)/hy  (3-49) 

where  wfi  = (w^3  + wg)/2  and  wA  = (^12  + w7)/2 

6w,x  = (6wd  - 6wc)/hx  (3-50) 

where  6wd  = (6w13  + 6w12)/2  and  6wc  = (6wg  + 6 w?)/2 
therefore 

w,y  6u,x  * <w8  ' w7  + w13  ‘ w12>  <{ul2 

- 6w?  + 6w13  - 6w8)/4hxhy  (3-51) 

thus  in  total  agreement  with  general  form  of  Equation  (3-47). 
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Mention  has  been  made  of  volume  referenced  to  a partic- 
ular derivative  approximation.  The  volume  is  defined  as  the 


area  h h times  the  ordinate  value  at  a given  grid/reference 
x y 

point.  By  treating  the  volume  as  a parallelopiped,  errors 
are  developed  in  the  total  integral  since  step  function 
expressions  are  created  with  respect  to  grid  arrangement. 

The  effects  are  more  noticeable  at  larger  mesh  size  (lower 
number  of  degrees  of  freedom) . 

The  above  paragraphs  have  shown  how  the  finite  differ- 
ence approximations  replace  the  partial  derivatives  within 
a virtual  work  expression,  now  the  steps  leading  to  the 
final  buckling  equation  are  discussed.  The  first  step  is 
to  factor  all  common  expressions  that  relate  to  specific 
6Wij  s.  Thus 


M N 

E E (C. . + A..  N ) 5w  . = 0 
i=l  j=l  13  13  x 13 


(3-52) 


where  is  the  only  externally  applied  force  considered. 


1 

. = £ £2  £ M - 2 £ M 

] sy . h . x.  x.. 

1 Jj  x l+l  l+l,  j 1 13 

1 

+ Cx  Mx  + 5X  IF  Cy  My 

xi-l  xi-l,jJ  xi  y L y3+1  i'i- 

" 2 My  + 5y  M + 

Yj  yij  yj-l  yi»j~lJ  hxhy 

nx  hv  M . ” '"'x  ^v.^xv. 

_xi-l  yj-l  xyi-l, j-1  xi-l  yj  xyi-l, j 

-nn  m +nnM 

Wi  V yj  xyij  J 


(3-53) 
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The  reader  is  directed  to  Appendix  B for  the  development  of 


several  terms  comprising  coefficients  C. . and  A. ..  It  be- 

i]  ID 

comes  obvious  that  in  Equation  (3-52) , the  only  way  a zero 
result  can  be  obtained  is  for 


+ A. 


i3 


N 

x_ 


0 


(3-55) 


since  6w^j  are  arbitrary  values  throughout  the  domain. 

There  is  a second  ordering  step  required  for  evaluat- 
ing a buckling  quantity.  In  Equation  (3-55)  the  displace- 
ment terms  (w_s)  occur  (as  can  be  observed  in  Appendix  B) 
indiscriminately.  It  therefore  becomes  necessary  to  factor 
and  order  all  the  w^j's  leading  to  an  eigenvector  expres- 
sion. 


(3-56) 


The  E . . matrix  contains  the  quantity  N (in-plane  load) 

1 J X 

which  is  the  eigenvalue.  Reference  [1]  discusses  the 
numerical  approach  used  in  the  final  reorientation  of  the 
displacements.  This  approach  is  called  the  marching  tech- 
nique and  is  developed  in  Appendix  C of  Reference  [1] . 


Finite  Difference  Approximation  of  the 
Equilibrium  Equation 

It  has  been  stated  in  Sec  II  Equation  (2-4)  that  the 
equation  for  buckling  is  a fourth  order  partial  diffential 
equation.  The  form  of  this  equation  is  restated  here  for 
convenience. 


The  finite  difference  approximation  to  this  equation  is 
established  by  making  use  of  full-station  approximations 
which  were  developed  earlier  in  this  section.  This  dif- 
ference technique  can  be  found  in  several  references.  An 
excellent  reference  which  presents  a detailed  example  of 
a plate  buckling  problem  is  "Structural  Analysis"  by 
Ghali  and  Neville  [26] . Thus  for  further  elaboration  on 
the  equilibrium  approach  the  reader  is  directed  to  the 
above  reference.  Yet,  for  continuity  the  finite  difference 
molecule  is  presented  in  Figure  3.7. 
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IV.  Numerical  Results 


Introduction 

Two  sets  of  boundary  conditions  were  investigated.  The 
following  terminology  is  applicable  to  both  the  simple  sup- 
port and  clamped  conditions. 

1.  Long  Plate  - This  method  employs  the  trigonometric 
finite  difference  approximations  by  choosing  the 
wavelength  parameter  X based  upon  the  buckle 
length  of  an  infinitely  long  plate  [1] . The  dif- 
ference approximations  are  then  substituted  into 
the  virtual  work  equation. 

2.  Conventional  - This  method  is  also  used  in  the 
virtual  work  equation  and  represents  the  polynomial 
approach  using  half-  and  full-station  finite  differ- 
ences. 


3. 


4. 


Equilibrium  - This  label  applies  to  the  full- 
station  polynomial  finite  difference  approximations 
of  the  equilibrium  equation. 

Lambda  - This  label  represents  the  values  used 
in  the  virtual  work  equation  when  the  buckle 
wavelength  is  not  based  on  the  buckle  length  of 
an  infinitely  long  plate 

Lamdba  - X^/A  (4.1) 

Values  of  Lambda  are  included  and  discussed  to 
present  a broader  appreciation  of  the  trigonometric 
finite  difference  capabilities. 


* 


Throughout  this  presentation  an  error  of  plus  or  minus 
five  percent  is  considered  acceptable  as  is  normally  the 
case  in  engineering  calculations  [27]. 

Various  Aspect  Ratios  - Simple  Support 

The  long  plate  theory  approach  gives  accurate  findings 
for  simply  supported  plates  under  in-plane  loading  if  the 
aspect  ratio  (A/B)  is  an  integer  (see  Fig  4.1).  The  reader 
should  notice, for  fractional  ratios,  accuracy  is  affected 
by  assumption  of  the  buckle  length  equal  to  the  width  of  the 
plate.  At  an  aspect  ratio  of  0.5,  the  plate  under  this 
assumption  is  more  flexible  while  for  a ratio  of  1.5,  the 
structure  is  stiffer.  Yet,  comparatively  speaking,  a long 
plate  assumption  incorporated  into  finite  difference  approxi- 
mations of  the  virtual  work  expression,  yields  results  more 
favorable  than  the  conventional  virtual  work  or  equilibrium 
equation  approach.  The  last  two  approaches  yield  about  the 
same  results.  One  should  also  notice  that  increasing  the 
degrees  of  freedom  (D.O.F.)  allows  a better  evaluation  of 
boundary  condition  effects  and  hence  less  error. 

Various  Aspect  Ratios  - Clamped 

Many  of  the  above  statements  become  appropo  for  clamped 
plates  (Fig  4.2).  Long  plate  theory  gives  the  best  results 
for  mid-range  aspect  ratios  while  conventional  virtual  work 
and  equilibrium  results  parallel  each  other.  It  appears  for 
aspect  ratios  1 to  4 at  least  81  D.O.F.  are  needed  to  keep 


accuracy  to  5 percent.  One  notices  at  an  aspect  ratio  of 
0.5  that  the  results  are  more  favorable  than  for  the  rest 
of  the  aspect  ratios.  For  this  aspect  ratio  we  have  a 
situation  in  which  the  boundary  conditions  are  more  realis 
tically  modeled  by  the  constant  mesh  arrangement  since  the 
plate  is  acting  similar  to  a wide  beam.  At  higher  aspect 
ratios  a constant  mesh  size  doesn't  properly  credit  the 
influences  of  the  boundary  effects.  Thus  to  obtain  5 per- 
cent accuracy  for  a buckling  clamped  plate  one  should  use 
at  least  81  D.O.F.  and  refine  the  mesh  size  near  the  bound 
aries. 

Degrees  of  Freedom  versus  Error  - Simple  Support 

Aspect  Ratio  0.5.  This  discussion  concerns  Fig  4.3 
and  the  rectangular  plate  in  this  case  is  similar  to  a 
beam  where  the  boundaries,  parallel  to  the  in-plane  load 
Nx»  are  far  away  from  the  load  itself.  In  assuming  long 
plate  theory  the  results  for  low  degrees  of  freedom  are 
very  good  in  fact  the  comparative  curves — equilibrium, 
conventional  and  Lambda  (>1.)  give  excellent  results.  In 
essence,  long  plate  theory  makes  use  of  the  fact  that  the 
buckled  length  is  equal  to  the  width  of  the  plate,  this 


can  be  illustrated  by  Fig  4. a 


Now  when  a value  of  1 is  selected  for  Lambda  (X^/A) , this 
assumes  a buckle  length  in  the  y-direction  equal  to  one- 
half  B (Fig  4.b).  This  comes  from  a basic  equation  (4.2) 
developed  in  Reference  [1]  to  calculate  the  parameters. 


A X 

i _5,  — D 

b • /T 


(4-2) 
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3 is  equal  to  the  aspect  ratio  for  simply  supported  plates 
and  1.5  when  the  plates  are  clamped  [28],  The  reason  the 
accuracy  isn't  as  good  as  long  plate  theory  is  that 
equation  (4-2)  for  Lambda=l  shows  the  two  buckle  lengths 
in  the  y-direction  rather  than  an  ideal  1 half-sine  wave 
for  the  length.  Considering  Lambda=0.5,  then  Ay/BT0.5=0.5 
yields  B=4Ay  thus  now  there  are  four  buckle  lengths  in  the 
y-direction  and  two  in  the  x-direction  hence  a stiff er 
plate.  This  is  realistically  pointing  out  what  Timoshenko 
[ 6 ] describes  graphically  (whenever  more  than  one  half-sine 
wave  is  assumed  in  the  y-direction  the  plate  is  stiffer) . 

The  curves  also  indicate  a band  of  Lambda  values  exist 
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from  one  to  infinity  (Conventional  method)  that  present 
excellent  results  at  low  D.O.F.  as  do  the  Equilibrium  and 
Long  Plate  methods,  while  for  values  of  Lambda  approaching 
zero,  the  results  become  erroneous  as  the  plate  becomes 
stiffer  (from  a buckling  viewpoint) . 

Aspect  Ratio  1.0.  Again  excellent  results  at  low 
D.O.F. , with  the  Long  Plate  method  showing  no  error  (Fig 
4.4)  as  was  expected  since  buckle  width  does  equal  buckle 
length.  A band  of  Lambda  values  greater  than  or  equal  to 
0.8  provides  good  results.  Aspect  ratios  2 and  3 also 
follow  these  trends  (Figs.  4.5  and  4.6). 

Aspect  Ratio  5.0.  In  this  instance,  the  curves  (Fig 
4.7)  have  inflection  points  at  lower  degrees  of  freedom. 

The  reader  should  remember,  especially  for  this  set  of 
curves,  the  concept  being  incorporated  into  the  analysis 
by  assuming  various  values  of  Lambda.  In  effect,  a func- 
tional value  for  the  displacement  is  assumed  at  buckling 
in  the  virtual  work  approach.  The  form  of  the  function 
selected  is  like  that  used  in  the  Rayleigh-Ritz  technique 
which,  in  essence,  requires  at  least  the  satisfaction  of 
the  kinematic  boundary  conditions  for  rapid  convergence. 

It  becomes  apparent  that  certain  Lambda  ratios  do  not 
provide  functions  that  meet  the  boundary  conditions  properly, 

(see  Fig  4.c),  thus  one  can  expect  greater  inaccuracy  at 
lower  D.O.F.  The  Long  Plate  method,  however,  gives 
excellent  results  at  low  D.O.F.  because  the  boundary 
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At  higher  D.O.F.  (64)  all  curves  present  monotonically 
decreasing  error  values.  This  is  expected  since  the 
assumed  function  becomes  secondary  with  its  inclusion  into 
the  virtual  work  equation. 

Degrees  of  Freedom  versus  Error  - Clamped 

Aspect  Ratio  0.5.  This  discussion  concerns  Fig  4.8. 
Again  we  have  the  conventional  virtual  work  and  equilibrium 
method  providing  almost  identical  results  while  the  long 
plate  theory  in  the  trigonometric  approximations  gives 


slightly  better  results.  It  is  also  noticed  that  a "good" 
selection  of  Lambda  (0.65)  produces  better  results  at  low 
D.O.F.  This  is  expected  and  explained  by  the  fundamental 
assumption  of  Long  Plate  Theory  where  the  width  equals  the 
buckle  length.  In  reality,  however,  the  clamped  boundary 
prohibits  a full  half-sine  wave  (see  Fig  4.e).  Thus,  with 
Lambda=0.65,  X^=.975b  which  is  more  realistic.  Aspect  1 
and  2 present  the  same  tendencies,  only  higher  D.O.F.  are 
needed  for  better  accuracy. 


Aspect  Ratio  3.0.  Concerning  the  interpretation  of 
the  Lambda  band,  the  relative  closeness  of  the  Conventional 
and  Equilibrium  methods,  and  the  ability  of  the  trigonometric 
finite  difference  method  employing  Long  Plate  theory  to  be 
slightly  better  than  the  Conventional  and  Equilibrium 
methods,  (Fig  4.11)  one  can  use  the  previously  stated 
observations.  Inflection  points  again  occur  at  low  degrees 
of  freedom  showing  a tendency  to  be  oscillatory.  An 
explanation  of  this  phenomenon  was  previously  presented. 

One  can  observe  that  even  with  the  conventional  approach, 
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potential  energy  wise,  its  more  accurate  to  select  25  D.O.F. 
versus  36.  This  peculiarity  is  attributed  to  the  failure 
of  the  selected  function  to  appropriately  consider  the 
boundary  conditions.  Salvadori  [29]  in  analyzing  equilibrium 
methods  experienced  oscillations  similar  to  these  for  certain 
cases.  It  should  be  noted  by  the  reader  that  when  inflec- 
tion points  do  occur,  long  plate  curves  are  not  effected. 

This  was  true  in  the  simple  support  case  (Fig  4.7)  and  also 
is  true  in  the  clamped  cases  (Fig  4.11  and  4.12).  The 
assumption  of  buckle  width  equal  to  the  width  of  the  plate 
seems  to  incorporate  boundary  conditions  more  accurately 
(like  the  simple  support  cases) . This  again  indicates 
that  the  inflection  point  is  probably  being  caused  by  poor 
approximation  of  boundary  conditions  at  low  D.O.F.  as  well 
as  an  insufficient  number  of  node  points  adjacent  to  the 
boundary  which  was  previously  pointed  out. 

Aspect  Ratio  5.0.  Oscillation  is  now  apparent  and 
troublesome  for  all  curves  (except  long  plate)  in  Fig  4.12. 
Higher  D.O.F.  are  needed  to  maintain  accuracy  and  this  is 
expected  since  again  we  have  increased  the  mesh  size  at 
low  D.O.F.  and  the  finite  difference  approximations  don't 
represent  the  boundaries  accurately.  Previous  Lambda 
interpretations  are  still  consistent  for  this  long  clamped 
plate  as  well  as  the  close  proximity  of  the  conventional 
and  equilibrium  curves.  The  value  of  the  trigonometric 
approach  can  be  appreciated  here  when  it  is  apparent  a 
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Lambda  value  of  0.10  would  have  given  accurate  results 
around  64  D.O.F.  compared  to  81  or  higher  for  all  other 
selections . 


FIGURE  4.2  FINITE  DIFFERENCE  APPROXIMATIONS 
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FIGURE  4.3  SIMPLE  SUPPORTED  PLATE  R/B=0.6 
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V.  Conclusions 

This  thesis  has  compared  results  using  expressions  of 
virtual  work  and  equilibrium,  employing  finite  difference 
approximations,  over  a broad  spectrum  of  aspect  ratios  for 
a buckling  plate  which  was  both  simply  supported  and  clamped. 
Little  difference  in  accuracy  was  found  between  results 
obtained  from  conventional  virtual  work  equation  and  the 
equilibrium  equation.  The  virtual  work  method  can  be 
improved  by  incorporating  into  the  finite  difference  expres- 
sions a trigonometric  function  using  a buckled  shape  for  a 
very  long  plate.  In  many  instances  certain  selections  of 
Lambda  values  lead  to  more  accurate  answers  with  the  trigo- 
nometric approach. 

Except  for  an  aspect  ratio  of  5.0,  all  methods  using 
simply  supported  plates  were  accurate  and  reliable  at  all 
degrees  of  freedom  (25  or  above).  At  an  aspect  ratio  5.0, 
inflections  of  the  error  function  were  encountered  and  some 
methods  oscillated  between  stiff  and  flexible.  The  reason 
apparently  is  that  large  constant  mesh  sizes  are  not  able 
to  properly  model  the  boundary  conditions  for  long  simply 
supported  plates. 

For  clamped  plates,  buckling  coefficient  inaccuracies 
began  to  increase  for  aspect  ratios  as  low  as  1.0  and 
oscillation  tendencies  were  encountered  for  aspect  ratios 
of  3.0  and  5.0.  Thus,  the  author  feels  that  additional 
reasons  for  the  inaccuracies  and  oscillations  must  be  put 
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forth.  They  may  be  thought  of  as: 

1.  The  inaccuracy  encountered  by  approximating  the 
boundary  conditions  in  the  basic  virtual  work 
equation  with  an  insufficient  number  of  node 
points  adjacent  to  the  boundary, 
and  2.  The  choice  of  a displacement  function  (referred 

to  as  Lambda)  which  doesn't  satisfy  the  kinematic 
boundary  conditions  properly.  (It  should  be 
noticed  that  the  Lambda  selection  is  very  similar 
to  choosing  a displacement  function  for  the 
Rayleigh-Ritz  technique  and  though  this  thesis 
hasn't  addressed  all  the  ramifications  of  this 
comparison,  several  inaccuracies  created  through 
the  use  of  various  Lambda  ratios  have  been  men- 
tioned. ) 

In  essence,  rapid  convergence  can  be  achieved  by  select- 
ing a Lambda  ratio  which  satisfies  the  boundary  conditions. 

It  should  also  be  noted  that  since  the  function  selected  is 
incorporated  into  the  finite  difference  approximation, 
accuracy  can  be  maintained  by  using  more  degrees  of 
freedom. 
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Appendix  A 

Development  of  the  Virtual 
Work  Equation 


The  strain  energy  of  a two-dimensional  plate  problem 
in  terms  of  strain  is 


U = f \ [e2+e2+2ue  e +^-  e2  ~ldV 

Jv  2 (l-u2)  Lx  y x y 2 xyj 


(Al) 


where  the  nonlinear  strain-displacement  relationships 

e =u  + ^-(w  ) 2 

x ,x  2 ,x 


V\y+  I<w,y>2 


(A2) 


e =u  +v  +w  w 
xy  ,y  / x r x ,y 

Equations  (A2)  are  substituted  into  equation  (Al)  and  then 
integrated  over  the  thickness  (t) 

0 - C C fit  7 7^  [lu,x4‘“.x»2l2+^.y4'“.y»2l2 


+2u  [u 

# 

x4(w,x>2lIv,y4 

<w,y)2] 

l-u 
+ “5“ 

lu,y+v,x+w/xw,y) 

dzdxdy 

(A3) 

u(x,y,z) 

= "zw,x 

v(x#y,z) 

= -«w#y 

(A4) 

w(x,y)  = 

w (no  dependence  on  z 

for 

small  deflections  in  thin 
plates) 
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Like  terms  in  w are  combined 


SU  = Uo  71^77  Iw.xx+uu,yy1Sw,xx+TI  -^77  [\yy 

+uw  ~1  6w  T w <$w  + — a - * r w3 

,xxj  ,yy  6 (1+u)  ,xy  ,xy  2(l-\>2)  L *xx 

+w  xw2vl  6w  x+ — "fc  » T w3v+w  vw2x  1 6w  v dxdy  (A8) 

#x  ry  j #x  2 (l-uz)  L-  #y  /X-J 
If  the  higher  order  terms  in  equation  (A8)  are  neglected, 
then  the  internal  virtual  work  of  a plate  during  buckling 
(specifically  derived  from  nonlinear  strain-displacement 
relationships,  but  representing  linear  strain  expressions) 
becomes 

5“  “ LL  (Mx4u.x*+MySw,yy+2Mxy6w,xy)dxay  (A9) 


Mx  * Dll",xx+D12w,yy 


My  “ D12",xx+D22w,yy 


(A10) 


M = 2D, -W 
xy  66  ,xy 


D22  = Et3/12 (1-u2) 


= Et  / 24 (1+u) 


(All) 


The  potential  energy  of  the  external  forces  is  completely 
developed  in  Appendix  D of  [16]  (Equations  are  shown 
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below  for  completeness) . In  their  development  it  is  pointed 


out  that  the  middle  surface  stresses  (N  N N ) are  assumed 

X/  y,  xy 

to  remain  unchanged  in  the  course  of  the  plate's  deflection. 
This  also  implies  that  no  middle  surface  stretching  occurs 
during  deflection,  thus  completely  eliminating  any  nonlinear 
large  strain  effects.  The  total  potential  energy  caused  by 
in-plane  external  forces  is 

V = i-  [ f In  w2  +N  w2  +2N  w w “jdydx  (A12) 

2 J0J0lx  *x  y *y  xy  *x 

The  virtual  work  is  then 

Sv  = 7 Do  [2V,x6",x+2V.ySw-y+2V'.ySw.x 

+2Nxy\x{w,y]dydx 
which  reduces  to 

SV  ’ Do  <NxW , xS w , x+NyW , y5w , x+Nxyu , x5w , y+Nxyw , y Sw , x’ dydx 

(A13) 

Equations  (A9)  and  (A13)  are  combined. 

This  gives  the  virtual  work  equation  used  in  this  thesis. 
rfa(MxS“.xx+My{“,yy+2MxySU,xy)dxdy  = (Nx“.x5u,x 


+V.y6”.y+N*yw.ySu.>‘+Vw.xSw,y)dydx  (A14) 


Appendix  B 
Factoring  Techniques 


This  particular  Appendix  is  designed  to  show  how 
numerical  equations  are  ordered  starting  from  the 


expression 


M N 

E E (C.  . + A.  .N  ) 6w.  . = 0 


i=l  j=l 


i]  ID  x 13 


(3-52) 


An  example  using  several  grid  points  is  presented  to  allow 
the  reader  a fuller  appreciation  of  the  numerical  volume 
integration  (Fig  3.4).  For  completeness.  Equation  (3-40)  is 
repeated  here  in  moment  form. 


6U  = hxhy  I=i  tllMnij/fix*D12Myij/f‘xfi2 

[S’'i+1,j-24“i3+iwi-1.l]+5xi5yj[I)12MXi/KJi:K 
+ D22Myij/fij]  [Swi,  j+l-2Swij+{ui,  j-l] 
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’ 6wij>/hx_ 


L Li  [5yjnx1Nx<wi+i.i-wii,(4wi*i.i 
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I 


Lt 


For  simplicity  let  M = M = Di2  = Dfifi  = 0 thus 

yij  yij  00 


Equation  (Bl)  with  (Be  becomes 


M N 
E E 
i=l  j=l 


S'yj  [>"*>3  [Swi*i,  j-2Swij+6wi-i.J 


f2 


+ [_cyjnxiNx,ui+i,-wij)  <Swi+i,:f  Swij,/hx 


= 0 


(B3) 


If  one  rearranges  and  combines  like  terms  the  following  is 
obtained 
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E E 


a4 

/u  ’ 


/h  H„  n„  N (W<i, 


p’v  11  v ' '•v  ‘v  “v  '"i  +1  -i 
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13 
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(B4 ) 


Consider  the  terms  i=2,  j=2  identified  as  1 2,3  2 


EE'  and 

corresponding  to  grid  point  7.  From  our  boundary  Equations 
(3-41)  thru  (3-44) 


i=2,j=2 


k2  2 y2  Z x2 


z 2 = j,  =1  n„  =1 


(B5) 


Equation  (B4)  becomes 
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Thus,  as  each  term  is  evaluated  the  virtual  displacements 

can  gain  coefficients  as  can  be  observed  in  Equation  (B6) 

for  Sw  and  6w0.  After  evaluating  all  the  terms,  the 

equation  is  reordered  to  gain  a final  coefficient  matrix  A 

and  a virtual  displacement  vector  6w^  where  the  coefficient 

matrix  A is  composed  of  C. . + A. . N . 
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